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A sea swell and surf model is implemented, tested and 
evaluated on a micro computer (HP-9 8 4 5B) . The two-dimensional 
model includes spectral wind-wave generation in open water, 
and shallow water wave transformation over irregular 
topography. The model predicts surf zone width, breaker 
lines and types of breakers. Using change in momentum flux 
of gravity waves (radiation stress) as forcing, the model 
predicts current velocities within the surf zone. The model 
is evaluated for the conditions over constant depths and 
uniform sloping beaches. The numerical results are checked 
against accepted theory and field observations. The model is 
found to overbuild wind generated sea heights for a 30 kts 
wind but to give expected wave heights for a 15 kts wind. 

The model results compare well with observed nearshore wave 
heights but give poor location of breaking waves. The 
model's nearshore current calculation is found unsuitable for 
the HP-9845B due to computational instability and time 
requirements. The sources of the model's problems are 
identified, and recommendations are made for future improve- 
ments of the model. 
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I . INTRODUCTION 



A. MOTIVATION 

Wave and surf forecasting is a primary concern of the 
U.S. Navy in executing amphibious operations. The wave 
environment affects the larger vessels in deep water and 
landing craft approaching the beach. As waves propagate from 
deep to shallow water, they undergo shoali j, refraction, and 
eventually breaking, while forming a surf zone. The height 
and type of breakers in the surf zone determines the difficulty 
that landing craft will have in reaching the beach. The 
breaking waves can also generate nearshore currents; these 
shear currents potentially can turn a craft and cause it to 
be broached and swamped by breaking waves. In addition to 
amphibious units, coastal structures and pierside ships are 
subject to forces by waves in shallow water regions. 

The forecasting of wave conditions can be divided into 
two regimes; the nearshore, including the surf zone, and an 
offshore area. There are a number of numerical models avail- 
able which predict waves, shear currents and surf conditions 
separately. Propagating wave energy from deep water to the 
high water line on the beach requires a combination of comple- 
mentary models. The objective of this thesis is to convert a 
state-of-the-art computer program for use on a Hewlett-Packard 
9845B mini computer, to test the model against field observations 
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and to revise where necessary. The program will be submitted 
as a prototype wave and surf model for shipboard use in the 
Fleet . 

B. BACKGROUND OF WAVE AND SURF PREDICTION 

The earliest models of wave generation and shallow water 
transformations were for monochromatic regular waves. Refrac- 
tion was based on the wave ray method (Arthur, 1949) , where 
ray paths were calculated for individual waves with slowly 
varying wave speeds along the wave front. This method is 
analogous to ray path prediction used in optics for light 
rays. With the development of computers and a better under- 
standing of the physical processes, it became possible to 
consider the entire wave spectrum in generation and shallow 
water transformation. Karlsson (1969) introduced a direction 
variable to study wave spectral transformation for straight 
and parallel contours. Collins (1972) included bottom friction 
and evaluated the validity of one dimensional models. Noda, 
et al. (1974) used a relaxation finite-difference scheme to 
solve the stationary wave spectral transformation in shallow 
water with the inclusion of bottom dissipation and effects of 
local wind generation. Shiau and Wang (1977) studied wave 
transformation in shallow water and compared the results with 
field observations. Wang and Yang (1981) applied their model 
with the inclusion of bottom friction. Recently, Chen and 
Wang (1983) improved the model by considering non-stationary 
waves and interactions with currents. 
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C. CURRENT STATUS OF MODELING TECHNIQUE 

The current state-of-the-art for numerical modeling can 
be characterized by two numerical methods, finite difference 
and finite elements. Both methods have been applied to wave 
and surf problems. A finite element model, such as by Wu and 
Liu (19 85) , economically handles localized irregular topography 
and has good properties for nonlinear convergence. However, 
finite element schemes are highly complex and require a sophis- 
ticated understanding of the model characteristics to choose 
the best grid system. Finite difference models (Noda, 1974; 
Shiau and Wang, 1977) use rectangular grids in general, and 
the numerical analysis has been widely studied. The numerical 
procedure of finite difference models is self-explained. There- 
fore, a finite difference scheme is less demanding on the 
user's numerical skills and from an operational point of view, 
is more practical. 

D. MODEL OVERVIEW 

The model described below was developed by Wang and Chen 
(1983). It is a finite difference model which considers a 
wind generated spectrum, spectral and swell wave transforma- 
tions over moderately irregular bottom topography, and 
dissipation by breaking waves and bottom friction. The model 
calculates a wind generated and deep water wave spectrum in 
the ocean and predicts wave and surf conditions on the beach. 
There is also an option of direct input of a predetermined 
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wave into a surf area. In the nearshore computations, th -re 
are options for computing surf conditions and/or nearshore 
circulation. The surf and nearshore circulation output are 
the basis for predicting amphibious landing conditions. Model 
input includes wind direction and speed, initial ocean swell 
and dominating direction, bottom type, and water depths. All 
the inputs are entered for the deep water calculation and 
again serve as appropriate data for the shallow water. 

E. OBJECTIVES 

The objective of this thesis is to compile the numerical 
model by Wang and Chen (1983) in a user friendly version and 
test it on the HP-9845B mini computer. The program should be 
self-explanatory with internal error checks and simple correc- 
tion steps. The model will be tested against observations on 
several beaches to establish accuracy and range of application. 
The physical and numerical limitations of the model will be 
documented, explained and where possible improved. It should 
be a useful tool for the Shipboard Aerographer Mates. 
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II. THEORY AND STRUCTURE OF THE MODEL 



A. THE SEA SWELL AND SURF PROGRAM (SSSP) 

The sea, swell and surf program (SSSP) was developed by 
Wang and Chen (1983) and is a combination of several comple- 
mentary models. The sea and swell for open sea follows the 
models by Shiau and Wang (1977) and Chen and Wang (1982) . The 
nearshore current calculation is basec on a model by Birkemeier 
and Dalrymple (1976). 

The version of the SSSP most suitable for use on the HP- 
98458 consists of a main program which controls nine subrou- 
tines. A flow chart of the program is given in Figure 1. The 
subroutine format of the SSSP allows for easier test and 
evaluation by component, modification or program correction, 
and reduced memory requirements on the mini computer. The 
main disadvantage of the subroutine format is some reduction 
in calculation speed. 

B. REFRACTION OF SEA SWELL 

The wave direction is determined using conservation of 
wavenumber for a steady state process (Phillips, 1977): 

3(k cos 0) 9(k sin 0) _ A /1x 

9y ” 9x 1 ; 

where (x,y) are the on-offshore and alongshore directions, 
and 0 is the counter clockwise angle from the X axis to the 
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Figure 1 . Sea Surf and Swell 



direction of swell propagation and k is the wavenumber. It 
is assumed that (1) is valid in both open water and within 
the surf zone. 

Noda (1972) developed an iterative scheme to solve (1) 
for the wave angle. The convention for grid location of the 
numerical model is shown in Figure 2. The differential equa- 
tion is central-differenced in the X-direction and forwarded- 
differenced in the Y-direction. The finite difference form of 
the equation usee by Noda (1972) with weighting on the forward 
difference term added to increase stability (after Abbot, 

1979) is used here: 



0 . . 

i/l 



sin 1 {j^ — [ (t ) (k sin 0) i+1 • 1 

i/ j ,J 



- (1 - 2x) (k sin 0) 



i+lf j 



( 2 ) 



+ (r) <k sin e) i+1<j+1 



^<<k cos 9> i(j+1 - (k cos 9) i, j-i 1 1 > 



where weighting factor of t = 0.25 is used. The initial guess 

for the wave angle is calculated by applying Snell's law for 

parallel bottom contours. (This gives an exact solution in 

the case of a plans sloping beach.) The initial guess field 

of 0 . .is iterated until 0. .is within a half percent of 
i/l i/l 

the previous iteration. 
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LONGSHORE 




Figure 2. Differencing scheme used in the model. The 

grid points are identified by subscripts with 
the first subscript indicating position from 
the beach with increasing numbers farther 
offshore. The second subscript indicates 
position alongshore with increasing numbers to 
the right when the beach is viewed from a 
ship off the coast. 
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The wave number k. .is the positive root of the linear 

i/D 

wave dispersion relationship 

a = (gktanh kd) (3) 

where d is the water depth and a is the angular frequency. 

In the original SSSP model by Wang and Chen (1983) , the wave- 
number was solved for in (3) using an iterative Newton's method. 
However, to speed up the model, the wavenumber is calculated 
here using a 6th degree polynomial fit of (3) as given by 
Hunt (1979) : 

k = — (S + (1 + 0.66667S + 0.35550S 2 

/ gd 

+ 0.16084S 3 + 0.06320S 4 + 0.02174S 5 

(4) 

+ 0.00654S 6 + 0.00171S 7 + 0.00039S 8 
+ 0.00011S 9 ) _1 ) " 1 



where : 

S = o 2 (d/g) 

This non-iterative approximation resulted in an observed 90 
percent savings in computation time on the HP-9845B with the 
same accuracy. 
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C. WAVE HEIGHT TRANSFORMATION 

The information of wave direction is used to calculate 
an initial guess for the energy field. The initial energy at 
each grid point is calculated by assuming conservation of 
energy over straight and parallel contours. 



3(E C cos 0) 

2 

3x 



0 



(5) 



Integrating (5) leads to the initial energy field: 



E = 



C cos 0 

E _2° o 

o C cos 0 

g 



(6) 



where the subscript (o) refers to conditions at the offshore 
boundary and is the group velocity. The initial guess 
energy field is used to calculate the wave height field. 

The calculation of wave height is based on the energy 
balance equation including refraction, dissipation and growth. 
Three assumptions are made to speed calculation and reduce 
model complexity (Chen and Wang, 1983). First, it is assumed 
that wave energy within a particular spectral band stays in 
that band, which allows linear superposition of wavelets. The 
second assumption is that each frequency component can be 
described by a single mean direction, 0. Third, it is assumed 
that there are no wave-current interaction. With these 
assumptions, the steady state energy flux equation for each 
spectral band can be written ( Longue t-Higgins , and Stewart, 
1960, 1962) . 



16 



p g |^(S(f) C g cos 6) + p g |^(S(f)C g sin 9) 

= - e d (f) + e g (f) 



( 7 ) 



where : 



S(f) 

P 

g 

e d ( f) 



wave energy density at spectral frequency f 
density of the fluid 
acceleration due to gravity 

energy dissipation through bottom friction 
energy generation due to wind stress 



The energy generation term, e , is based on the mechanisms 

s 

used for wind energy transferred to the wave field described 
by Phillips (1957) and Miles ( 19 57 ; 1959a ,b ; 1962 ) . The Phillips 
mechanism is due to random atmospheric pressure fluctuations 
advected by wind resonantly interacting with a perturbed sea 
surface. The wave energy generated by the Phillips wind 
mechanism is 



where 



and 



S p = P ga 



,2,3 
4 7T k 

a = 2 3" P(k,a) 

P 9 



( 8 ) 



(9) 



a, = a/g 



( 10 ) 
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Barnett (1968) fitted the following empirical relationship 
for P based on field measurements: 



P(k,a) 



6.13 x 10 4 w 6 

2 2 
7T a 



V, 



v 2 + (k sin 6) 2 



3 



where ; 



v. 



v 2 + (k cos 6 - A) 2 



3 



( 11 ) 



6 



w 



v 

V 



A 

1 

2 



angle between wind and wave 
wind speed (meters/second) 
a/w 



0.3 



1.28 



0.52 



0.95 



This mechanism is important for wave initiation from calm 
water and results in a linear increase in energy with time 
(Hasselman, 1960). 

The Miles energy transfer mechanism is important only 
after waves have formed (LaBlond and Mysak, 1978). When wind 
blows over an irregular sea surface, flow separation can 
occur in the lee of the wave crest creating pressure differ- 
ences. Momentum is transferred to the waves as a result of 
these pressure forces. The instability mechanism of Miles' 
results in exponential wave growth (Barnett (1968) ; Hasselmann 
(1960) ) . 
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pgB [S ( f ) ] 



( 12 ) 



where : 



ra 



6 



As a ,W cos 6 
2 IT 1 c 



8 ] 



(13) 



s = ratio of air to water density 
A = 5 

B = .9 



The wind wave generation is then the sum of 



<V 



and 



(s m ) 



The wave growth is limited to a fully risen sea at each 
spectral band. The limits on development are based on Phillips 
(1957) equilibrium spectrum modified by Kitaigorodskii et al^ 
(1975) and Thornton (1977) for the shallow water condition: 

S(f,5) = B g 2 f" 5 H(6) T(W n ) (15) 



where : 



H ( 6 ) = 



(£) (cos 4 (5) ) 

J 7T 



(16) 



accounts for angular spreading, 



_ 9 2W I (W ) -1 

r(w ) = £ (w ) (l + — 

sinh[2W Z(W ) 
n n 



-} 



(17) 
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and 



W 



n 




(18) 



where 



nw n ) 



satisfies the condition 

Z Tanh (W Z) = 1 

n 



F (W n ) = 1 in deep water, and for W n less than 1, 

r<w„> = 



(19) 



( 20 ) 



and 



S (f , 6 ) = j qg d f -3 H ( 6 ) (21) 

_ o 

where q is a coefficient on the order of 10 ; a value of 

.073 is used in Wang's model. 

The bottom dissipation term, e^, is due to work done on 
the bottom by the bed shear stress. The bed shear stress is 
proportional to velocity squared and the rate of energy 
dissipation is proportional to the cube of the velocity. 
Hasselman and Collins (1968) proposed a quasi-linearized 
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expression for a frequency-de '.dent dissipation function: 



e d (f) = P c f |U b | u b (f) (22) 

where is a bed shear stress coefficient, u b (f) is the wave 
induced bottom velocity at a particular spectral component, 
and U b is the bottom velocity for the total flow field. 
Applying linear wave theory transfer functions relation the 
energy spectrum of the surface elevation to the velocity 
field, (22) can be expressed as: 

2 2 

e.(f) = pc f K — 2-^ ]S(f)U. (23) 

d f a 2 cosh 2 (kd) b 

where the term in the brackets is the spectral transfer 
function. The total bottom velocity is obtained by applying 
the spectral transfer function to the entire spectrum: 

(U ) = [ / -9 k S(£) df] (24) 

0 a cosn kd 

D. THE SURF MODEL 

This model describes the surf zone width, number of lines 
of breakers, breaker type, breaker angle, significant and 
maximum breaker height, and the effective surf. The calcula- 
tions are based on empirical formula given in the Shore Pro- 
tection Manual (U.S. Corps of Engineers, 1977). The input 
wave is monochromatic, and requires only significant wave 
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height and period. Initial wave height and direction are 
determined for the entire model grid using swell direction 
and wave height calculations described in Chapter II, sections 
B and C. The depth of the first breaker is determined by 
comparing the calculated wave height with the incipient break- 
ing criteria (Noda et ad., 1974; Thornton and Guza, 1983) . 



The width of the surf zone, X^, is then calculated using 
the relationship between breaker height, breaker depth, wave 
period, and bottom slope (Weggel, 1972) . 



H 



_b 

L 



0.43 Tanh (— ) 
Xj 



(25) 




b(m) - a(m) 

T 




(26) 



assuming 




c + dX, 



b 



(26a) 




e + fX. 



b 



(26b) 



Solve : 



(a'bf)X^ + (a'ed + a ' cf - b'f + d)X^ + (a'ce + c-b'e) 



b 



0 



(26c) 
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The coefficients c, d, e and f are determined by using 
standard slope intercept techniques. The positive root of 
this quadratic equation is the width of the surf zone. 



a(m) = 1.36(1 - e ^ m ) 



b (m) 



1.56 



1 + e 



-19.5m 



a' = 



a (m) 

~IT~ 



b' = b(m) 



(27) 



(28) 



where : 



m = bottom slope between grid points 

The first breaking wave height, denoted as HBK in the 
program, is then determined using the surf zone width by 
interpolating between grids. 

H I J " H I-1 J 

HBK = ((H ((2-1) + (H H ) (I-D) + ( lf ■ Ax - ~ ) x 

(29) 

where H is the original wave height from energy considerations 
only and I and 1-1 are the grids that bracket the deepest 
breaker in the on-offshore direction. 

The depth of the first breaker, DBK, is calculated also 
using the surf zone width an -he interpolation between grids: 
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DBK 



- "“i.j" 2 - 1 ) + + ( Dl,J -T x~ 1 ,J 



)X 

( 30 ) 



where D is the bottom depth in the program. 

Additional lines of breakers shoreward of the initial 
breaker are calculated by: 



H b next 



b next 



1 + a' d 



b next 



(31) 



where 



b next 



h I,J + 



(32) 



The angle of the breaking wave (Angz) is calculated by 
interpolating the distance between grid points and the linear 
change in the wave direction: 



Ang - Ang 

Angz = Ang + ( • ' ~ K — ) (X - ( 1-1) Ax) (33) 



where Ang is the calculated wave direction at the grid point. 
The type of breaker is based on the surf parameter: 



m 



(H /L ) 
o o 



T/2 



(34) 



where : 
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m = bottom slope 

H = deep water wave height 
2 

L q = 1.56T , deep water wave length 

Breaker type is specified by the following criteria: 

> 2 surging or collapsing 

0.4 < < 2 plunging 

< 0.4 spilling 

Effective surf is calculated in accordance with 
COMNAVSURFPAC/COMNAVSURFLANT Instruction 3840.1 for use in 
naval amphibious operations. "Effective surf height is a 
numerical value, expressed in feet, and provides a guide for 
judging the feasibility of landing operations for any given 
landing craft under existing or forecast conditions." Effec- 
tive Surf is a function of: 

1. Breaker Height 

2. Breaker Period 

3. Breaker Type 

4. Littoral Current 

5. Wind direction and velocity 

6 . Secondary wave height 

7. Breaker angle 

The parameters for the empirical calculation of effective 
surf are calculated in the surf condition part of the model 
except for littoral current. The nearshore current circulation 
in the model is described in the next section. 
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E. NEARSHORE CIRCULATION MODEL 



The nearshore circulation is generated by changes in momen- 
tum due to breaking waves within the surf zone and by wind 
induced surface stress. The nearshore circulation model is 
based on the depth integrated and time averaged equations of 
motion by Phillips (1977) . The breaking wave-induced currents 
are calculated by alternately solving the momentum and the 
continuity equations. The equations are solved numerically 
with the finite differences scheme taken directly from Kirby 
and Dalrymple (1982) which is an updated version of Birkermeier 
and Dalrymple (1977) . The x and y momentum equations are 
integrated over the depth and time averaged over one or several 
wave periods. The equations contain bottom shear stress, mean 
surface shear stress and excess mean momentum stress due to 
wave action. To simplify the model, the formulation neglects 
the nonlinear mean momentum flux terms and the lateral (hori- 
zontal) transfer of turbulent momentum. The continuity equation 
is given by: 



where rj" is the free surface elevation, d is the still water 
depth, (u,v) are the (x,y) horizontal velocity components and 
the overbar indicates time averaging. The horizontal momentum 
equations are 



§t + L tu(d + n>] + Livid + n)J 



0 



(35) 



3u 

3t 




(36) 
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1 



(37) 



3v 

at 



g I 11 - 
9 3y 



as 



■( 



yx 



p ( d+n) 



— . 3x 



as 

_ n. . 

ay 



sy 



T by> 



The S , S , S , S are the mean momentum flux due to the 
xx xy yx ' yy 

wave motion (commonly referred to as radiation stresses) . 



T , T 


are 


the 


surface 


shear stresses due to wind 


sx' 


sy 








T, , t, are 

bx' by 


the 


bottom 


stresses . 



The radiation stresses are given by (Longuet-Higgins and 
Stewart, 1964) : 



S = E[(2n - l/2)cos^9 + (n - l/2)sin^9] 

XX 



(38) 



S yy = E[(2n - l/2)sin^0 + (n - l/2)cos^0] (39) 



S = S 
xy yx 



2 n sin( 20 ) 



(40) 



1 2 

where E is wave energy = 3 -pgH ; n = C /C. 

« g 

The surface wind shear stresses are calculated (Van Dorn, 
1953) by: 



x sx = P K|W|W X (41) 

x sy = P K|w|W y (42) 

where W , W are the wind speed in the (x,y) directions, and 
x y 

_ g 

K is the wind stress coefficient of order 
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The bottom shear stresses are linearized by assuming weak 



mean currents (U,V) (LaBlond and Mysak, 1978) . 

4 rT 

x, = p — c c u U 

bx tt f m 

2 

T. = P — Cj. U V 

by it f m 



(43) 



(44) 



where is the Darcy-Weisbach bottom friction factor and 
u^ is the maximum wave orbital velocity at the bed: 



H g 

u m 2 sinh kd 



(45) 



The maximum time step At is based on the linear stability 
criterion (Kirby and Dalrymple , 1982) : 

At < \/((Ax) 2 + (Ay) 2 /2gd) (46) 



However, in practice, Kirby and Dalrymple (1982) recommend 
using a quarter of this value. The actual At used in the 
program is 



At 



1 (Ax 
3 



+ Ay 2 ) 

172 “ 



(gd) 



1/2 



(47) 



The finite difference form of (36) and (37) are given 

by : 
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where : 



d = 



^ T T — 1 ^T_1 T— 1 

( — ■= ) is the average depth in 

the X-direction 
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where : 
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The finite difference form of (35) solving for n is 
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I,J 
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The new water depth is then given by: 



I,J 



old _new 
(d T _ - n, ,) + n. 



I,J 'I,J 



I,J 



The nearshore circulation solutions are "spun up" to 
steady-state condition. The model is started from state of 
rest with zero wave field present (Kirby and Dalrymple, 1982) . 
To reduce the effect of seiching or instability, the vave 
height H at the offshore grid is gradually brought up to 
its full value using: 



H 



H q tanh (— ) 



(51) 



where : 

t = model time 
T g = arbitrary fixed period 
H q = full height of boundary wave 
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The hyperbolic tangent function is used to control the 
growth of the boundary wave during an iterative process. The 
quantity [2t/T ] increases linearly while the tanh function 

W 

increases smoothly toward unity. 
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III. MODEL APPLICATION ON THE HEWLETT-PACKARD 9845B 



A. CONVERSION INTO BASIC 

The U.S. Navy contracted for the development of the wave 
and surf model which was written in FORTRAN. In its original 
form, the model was written for a large non-portable computer. 

To make the model useful tactically for mobile U.S. Navy forces, 
the program is rewritten specifically for deployed micro 
computers like the HP-9 84 5B. It is required that the model 
be written in a user friendly format with common naval conven- 
tions and terms. Error checking and input review had to be 
programmed into the model to be operationally acceptable to 
the U.S. Navy. The Wave and Surf Model had to run on the 
HP-9845B to give reliable, useful and easily understood output. 

B. MODIFICATIONS IN THE WAVE AND SURF MODEL FOR THE HP-9845B 

The original contractor delivered a FORTRAN version of 

the wave and surf model written in a format of "stand alone 
subroutines." This version was suitable for conversion to 
the HP-9845B because its size was within the micro computer's 
memory requirements. The subroutine format also allowed for 
easier testing of the model. The main change made in the 
program involved the use of common statements. The FORTRAN 
version made extensive use of labeled common blocks. The HP 
basic only allows one common statement. This minor problem 
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was solved by increasing the size of variable? n the call 
statement to most subroutines. Only the most ^iten used 
general variables were used in the single common statement. 

The language of the input statements is changed to agree 
with naval conventions. Wind and swell are described by the 
direction they arrive from. The beach is described from the 
point of view of a ship offshore. 

The interactive subroutines are expanded with lines of 
code that: 

1) explain the input needed 

2) ask for the data input 

3) error check the data input before accepting or 
rejecting and requesting correct input. 

The program was expanded in the calculation code with 
explicitly assigned variables for counters and matrices. For 
example, effective surf was stored into a unique array rather 
than sharing an array used for several other calculations. 
Empirical values from the COMNAVSURFPAC/COMNAVSURFLANT Instruc- 
tion 3840.1 were put into an array rather than use the FORTRAN 
version data fitting calculation. The additional memory 
space for these arrays was traded for explicitly assigned 
data and clearer programming. 

Labels were placed before every significant calculation 
in the model code. Spaces were added between code to separate 
functional groups. Indentations of the code were made to show 
loops and double loops. These changes added to the readability 
of the program and made the model logic easier to follow. 
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In addition, changes/corrections were made to the sea 
surf and swell model for undefined constants, incorrectly 
dimensioned arrays, logic errors and improperly defined 
angles. Some of these changes were made with the cooperation 
of the contractor. In the program code, each change is 
documented . 
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IV. DISCUSSION 



A. USABILITY 

The sea, surf and swell model is programed to be inter- 
active in language a sailor can understand. The units for 
wave height and bottom depth entires can be either meters or 
feet. Grid spacing can be in nautical miles only for the 
open water version. However, the nearshore model allows 
either feet or meters. All entries are requested with a 
format example. Where possible, the format example is in 
standard World Meteorological Organization code, which should 
be familiar to all U.S. Navy Aerographers Mates. Suggested 
input values, or default values, are offered where appropriate 
to make the model more "user friendly." 

B. DATA INPUT REQUIREMENTS 

The quality of the input data is the most important factor 
determining accurate wave and surf forecasts. The bathymetry 
begins to affect sea and swell when the depth is less than half 
a wave length. For the same depth of water, the effect of 
the bottom will be strongest on longer period waves. The 
more shallow the water, the greater the effect will be on all 
waves. Accurate bathymetric data is needed in deeper sections 
of the model grid, but becomes increasingly important for 
shallow waters. The wave period and depth is used to compute 
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wavelength and phase speed. These parameters determine 
forecast wave direction, energy, surf type and currents 
generated in the surf zone. For the same reason, initial 
wave direction is as important as the wind direction . The 
bottom type input influences energy dissipation and surf zone 
currents. The surf model calculation is least sensitive to 
the bottom type. 

The most difficult, and clearly the most important input, 
is the offshore wave input. Careful consideration must be 
given to the many sources of deep water waves. The open 
model sums the energy from all deep water wave inputs. The 
wave spectrum input must be accurately described. The ocean 
can be nearly flat concealing long period swell that can 
generate catastrophic surf on beaches. The ocean can have 
fully developed seas generated by continuous wind. In a 
fully developed sea, the energy of the waves is distributed 
over a wide range of periods that must be entered into the 
model for accurate surf forecasts. Only with a sound knowledge 
of deep water wave conditions can the model be used to fore- 
cast shallow water conditions. 

The open water part of the model has the most options for 
wave energy input, which include: 

1. A single narrow band swell. 

2. A significant wave height and period which can be 
expanded in terms of a Bretschneider spectrum. 

3 . Specific energy levels can be entered at appropriate 
frequencies such as from the Navy Spectral Ocean Wave 
Model . 
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It is important to note that the output from this model is 
a single significant wave height and period moving in the 
direction of maximum energy. This could be a poor represen- 
tation if the seas were large and from diverse directions. 

The surf and nearshore circulation parts of the model 
allow only a monochromatic wave input. This input could be 
directly from the open water calculation or an input of the 
user's choice. If there is more than one dominant wave group 
approaching the beach, the model representation of the surf 
zone will be incomplete. 

C. SELECTING THE GRID SIZE 

The grid size and the distance between grid points will 
determine the model's ability to resolve the wave climate. 

The open water model has grid distances in nautical miles. 

The grid mesh should cover the area from deep water waves 
into shallow water, outside the surf zone. The grid size for 
the surf and nearshore is measured in feet, and the selection 
is critical to resolving the surf zone. Kirby and Dalrymple 
(1982) outline the following requirements for the nearshore 
circulation grid. 

1. The grid must extend offshore far enough to remove the 
offshore region of the domain from the influence of 
currents driven by the surf zone, and to allow for 
the specification of a uniform longshore depth which 
will not significantly alter the wave refraction 
results in the nearshore. 

2. The grid mesh must be fine enough to resolve the surf 
zone adequately. 

3. In the event that the effect of a single physical 
feature isolated in the longshore direction is to be 
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modeled, the longshore extent of the grid must be 
sufficiently large to isolate the physical system 
from the effect of images created by the longshore 
periodicity requirement. 

Experience running the model suggests that a minimum of 
3 grid points must be inside the surf zone. This would not 
be possible for waves breaking directly on the beach. A 
further constraint is the number of grids in this model (10 
offshore rows and 13 columns alongshore) due to the limita- 
tions of speed and storage area of the micro computer. 
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V. MODEL TEST AND EVALUATION 



A. DEEP WATER WIND GENERATED WAVES COMPARED WITH JONSWAP 

CURVES 

The model approximates the area under the spectral density 
curve (wave energy) with rectangles. The more rectangles 
used, the better the approximation will be, but computational 
time also increases. The open water model requires approxi- 
mately 45 minutes to calculate each frequency band on the 
HP-9 84 5B. Narrow band swell waves might be closely approximated 
with a single energy band and a band width of .01 Hz. A 
broad band wave spectrum, such as wind waves, needs to be 
represented by multiple energy bands. 

To test the wind generation model and the effects of 
bottom friction and band width selection, the model was run 
for constant depths of 100, 10, and 5 meters (Figs. 3, 4, and 
5) . The initial wave height is zero. The bottom friction 
coefficient chosen is .01 for a fine sand bottom. The grid 
spacing is 10 nautical miles resulting in a maximum fetch of 
90 nautical miles. The results are compared with the JONSWAP 
significant wave height as a function of fetch from Bishop 
(1983). The model results favorably agree with the JONSWAP 
curve considering the relative coarseness of the energy 
bands. For the 100 m depths, bed friction was insignificant 
in changing the wave height as expected. 
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Figure 3. Significant wave height generated by a 15 kt 
wind in 100 m depth. Model computed values 
(solid line) are the sum of spectral contri- 
butions at .05, .10, .15 and .20 Hz frequency 

bands. The dashed curve represents wave 
height generated by the JONSWAP spectrum as a 
function of fetch and wind curves in deep water. 
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Fig. 4. Significant wave height generated by a 15 kt 

wing in 10 meter depth. Model computed values 
(solid line) are the sum of spectral contribu- 
tions at .05, .10, .15 and .20 Hz frequency 
bands. The dashed curve represents wave height 
generated by JONSWAP spectra as a function of 
fetch and wind in deep water. 
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Figure 5. Significant wave height generated by a 15 kt 
wind in 5 meter depth. Model computed values 
(solid line) are the sum of spectral contribu- 
tions at .05, .10, .15 and .20 Hz frequency 

bands. The dashed curve represents wave 
height generated by JONSWAP spectra as a 
function of fetch and wind in deep water. 
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The model was run for the same input with the flat bottom 
depth changed to 10 meters. The effect of friction only 
slightly reduced the fully arisen significant wave heights. 
However, the model generated fetch limited wave heights were 
approximately half a foot higher than the JONSWAP results. 

The increase in the height of the "shallow water waves" was 
probably a result of the change to a shallow water saturation 
spectrum. In shallow water the slope of the energy saturation 
spectrum with respect to frequency is -3 (Thornton, 1977) . 

In deep water, the slope is -5 (Phillips, 1966) . This results 
in relatively more wave energy at higher frequencies for the 
shallow water spectra. 

To test the bottom friction in the model, a run was made 
with the same input as before except the depth was changed 
to 5 meters. The significant wave height curve shows a slow 
decrease after the first grid point. The model calculates 
larger frictional dissipation as the down grid distance in- 
creases. The model nearly reached a constant wave height, 
balancing generation and dissipation forces after 80 nautical 
miles. The values of predicted wave heights by the JONSWAP 
spectrum agree with the present model within 10 nautical miles, 
but tend to be larger for a longer fetch. For an open sea in 
a shallow region, the JONSWAP spectrum ignoring bottom friction 
effect gives a wave height four times larger than the pre- 
dicted by the present model; this result emphasizes the 
importance of including bottom friction in shallow water wave 
models . 
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The wind generation part of the model was again tested 
for 30 kt winds blowing over the same area. The significant 
wave heights computed by the model increase rapidly as the 
downwind fetch is longer. The model predicts a near fully 
arisen sea after only 50 nautical miles. The JONSWAP model 
predicts approximately linearly increasing wave height values 
for 30 kt winds until reaching fully arisen seas after 200 
nautical miles of fetch. The large difference between the 
model calculations at 30 kts and the JONSWAP curves suggest 
the model over-builds seas for this wind speed and grid spacing. 

The energy transport calculation in the model was found 
dependent on grid spacing and is invalid for large grid 
spaces and higher frequencies. Negative energies resulted for 
input at frequencies greater than .3 Hz with a grid spacing 
of 10 nautical miles. By decreasing Ax = Ay to one nautical 
mile, the model can compute energy transport for .4 Hz and 
higher. A "fix up" in the program is made by setting any 
negative energy values to zero. The inability of the model 
to compute energy transformation for the short waves should 
be noted, even though this is the low energy side of the 
spectra . 

B. NEARSHORE CIRCULATION ON A SLOPING BEACH 

The nearshore current prediction part of the model was 
tested on a uniform sloping beach. Since the depth is uni- 
form in the alongshore (Y) direction, the wave height, wave 
direction, current, wave set-up and surf should also be 
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* Figure 6. Significant wave height generated by a 30 kt 

wind in 100 meter water depth. Model computed 
.values (solid line) are the sum of spectral 
contributions at .05, .06, .07, .08, .09, .10, 
.11, .12, .13 and .14 Hz. The dashed curve 

represents wave height generated by JONSWAP 
spectra as a function of fetch and wind in 
deep water. 
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uniform in the Y direction. The calculated wave directions 
were constant in the Y direction, the wave height and surf 
changed a few percent but the variations of current and wave 
set-up in the Y direction deviate more than 20 percent. 

The wave height and surf were slightly inconsistent with ex- 
pected results and this is thought to be caused by the finite 
differencing scheme producing errors from the boundaries. The 
larger variation in the current and wave set-up values is 
probably due to the inconsistent mass flux calculations. The 
boundary errors induced in the wave calculations resulted in 
a nonuniform forcing term, and, therefore, they may not be 
able to describe a constant current. The model is expected to 
give best results in the interior of the computational domain 
for wave and surf prediction. 

The nearshore current calculation portion of the model 
was compared with data acquired at Santa Barbara on 4 February 
1980. Initially the currents calculated were one order of 
magnitude smaller than the observations. 

Consequently, three basic problems with the surf model 
were identified: 

1. The model computed only 36 seconds of surf after 
beginning with still water. The contractor's program 
calculated only 10 time steps, or 3.6 seconds. This 
length of time was too short to even approach a 
steady state current. Kirby and Dalrymple (1982) 
always ran their nearly identical version of this 
model for a 1000 time steps representing more than 

5 minutes real time. 

2. The resolution was too coarse to describe the current 
field. The finite difference scheme is unable to 
compute currents nearer than 2 gird spaces from a 
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depth of zero (at the beach presumably) . The surf 
zone for this test was closer to the beach than 2Ax 
or 20 meters. 

3. The changes in radiation stress, which forces the 

current, were very small across the breaker line where 
it should be the largest. The current would remain 
too small even with items 1 and 2 corrected due to 
weak forcing. Realistic radiation stress was not 
being calculated because the model gave incorrect 
wave heights for breaking in the surf zone. This 
correction was made. 

The model was run again for comparison with Santa Barbara 
data. The grid size was reduced to Ax = 2.5 meters and 
Ay = 5 m. The time step for this small grid was approximately 
.6 seconds. The HP-9845B was allowed to calculate for 381 
time steps which took 9 hours. This calculation represented 
5.28 minutes real time for waves entering still water. The 
maximum stable time step is a function of the grid size 



TMAX = 1/3 

gd 



The current was approaching a steady state with a very slow 
rate of convergence (Fig. 7) . It would take more than 23 
hours of CPU time on the HP-9 84 5B to run the model for a 
1000 steps as done by Kirby and Dalrymple (1983). Attempts 
to increase the time step and decrease the CPU time resulted 
in numerical instability (Figure 8) . The model was run using 
the same input except the At was increased 5 percent each 
time step until reaching a maximum of one second. Instability 
began to occur after 20 seconds with values oscillating around 
the stable solutions. After 25 seconds the unstable model 
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Figure 7. Current generated by the surf as predicted 
by the model as a function of "spin up" 
time. A total of 381 time steps were made 
representing 316 seconds of forcing by the 
surf. The dashed line is an extrapolation 
of the model's calculation. 
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Figure 8. Surf zone current modeled for 4 Feb 1980 

off Santa Barbara, California. The stable 
increase in current is shown with the 
solid line. The unstable calculation is 
shown with the heavy dashed line. 
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of the current became large and changed sign. The surf model 
is very slow to converge to a solution and unstable for even 
slightly longer time steps. 

C. SURF CONDITION 

The surf zone wave heights and directions calculated by 
the model were compared with observations. The comparisons 
were made with observations for two California beaches at 
Santa Barbara and Torrey Pines. These beaches are nearly 
plain sloping. The offshore boundary condition of wave 
height, period and direction are given in Table I. The 
beach material for the model friction coefficient selection 
was fine sand (Cf = .01) . The modeled and observed wave 
direction agreed, within round-off, to the nearest degree 
for all days. The waves turned toward normal while crossing 
the surf zone exactly as expected. 

Model comparisons for observations on 3 February are shown 
in Figure (9) which shows the model slightly under-predicts 
the wave heights but follows the general change across the 
surf zone. The breaker line was forecast to be 14.9 meters 
from the shore with a height of .68 meters. The observed 
breakers began more than 40 meters from the beach with a 
height of .68 meters. The model incorrectly forecast one 
line of breakers when at least two were observed. 

The 4 February results (Figure (10)) show reasonable 
prediction of wave height. But using the wave height 
information, the model did not predict the location or height 
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TABLE I 







Offshore 


Boundary Conditions 




Date 




Depth 


H 

rms 


T Wave 


direction 






(m) 


(m) 


(sec) 


(°) 


Santa 


Barbara, Leadbetter 


Beach, California 




3 Feb 


80 


4.1 


.55 


14.3 


7.8 


4 Feb 


80 


4.2 


.56 


14.3 


9.0 


5 Feb 


80 


4.1 


.45 


12.8 


8.4 


Torrey 


■ Pines Beach, 


California 




4 Nov 


78 


6.7 


. 35 


14.3 


0 


10 Nov 


78 


6.6 


.68 


15.9 


0 
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Figure 9. Santa Barbara surf zone wave height on 3 
February 1980. Model computed values are 
represented by the solid line, observed 
heights are dashed. 
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Figure 10. Santa Barbara surf zone wave height on 4 
February 1980. Model computed values are 
represented by the solid line, observed 
heights are dashed. 
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or breaker. The observation indicated two lines of breakers, 
the first near 60 meters and the second near 38 meters. The 
maximum measured wave height was .67 meters occurring 60 
meters from the beach. The model forecast only one line of 
breakers .84 meters high and 16 meters from the beach. 

The 5 February results (Fig. (11)) show the model almost 
forecast the surf height exactly at grid points. The model 
forecast one line of breakers 11 meters from the shore with 
a height of .61 meters. The observed breakers began about 
39 meters from shore with a height of nearly .60 meters. 

The model was also compared with surf zone observations 
from Torrey Pines State Beach, California for 4 and 10 
November 1978 (Table I) . The average angle was near normal 
to the shore. Fig. (12) shows the model input boundary wave 
on 4 November height was too small by about .05 meters (2 
inches) . This small error was due to interpolation between 
sensor locations. The model followed closely the observed 
wave heights between 160 and 240 meters, but it over-built 
the wave height inside 160 meters. The model placed the line 
of breakers 82 meters from shore with a height of .76 meters. 
Observed breakers were near 160 meters with a height of .48 
meters . 

The 10 November 1978 model of Torrey Pines was based on 
boundary conditions entered at 6.6 meters. Figure (13) shows 
uniform over-prediction of the wave height at grid points. 

The model forecast two lines of breakers, the first line 133 
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Figure 11. Santa Barbara surf zone wave height on 5 
February 1980. Model computed values are 
represented by the solid line, observed 
heights are dashed. 
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Figure 12. Torrey Pines surf zone wave height on 4 

November 1979 . Model computed values are 
represented by the solid line, observed 
heights are dashed. 



\ 
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Figure 13. Torrey Pines surf zone wave height on 10 
November 1979. Model computed values are 
represented by the solid line, observed 
heights are dashed. 
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meters from the beach with a height of 1.7 meters and the 
second 83 meters from the beach with a height of .75 meters. 
However, the observations showed no corresponding lines of 
breakers. The observed heights of the lines of breakers are 
.72 meters or less. 

The wave height calculation by the model is, in general, 
close to the observed wave height. The surf zone width 
based on (26) is incorrect in every test run. The resulting 
calculations of breaker position, angle and height are 
dependent on the width of the surf zone and therefore these 
forecasts are incorrect. 
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VI. SUMMARY AND CONCLUSIONS 



A. SUMMARY 

A numerical model describing surface gravity wave trans- 
formation in shallow water region is written in BASIC for use 
on the HP-9845B. The original version of this model was 
based on the work by Wang and Chen (19 83) . The model consists 
of an open water wave program and a surf program. The open 
water model considers a wave energy spectrum generated by 
winds; the surf model predicts a shallow water wave field and 
the shear currents induced by wave breaking on an arbitrary 
bottom topography. The model is tested for cases of constant 
depth, a uniformly sloping beach and actual beaches. The 
application of the model to a field beach indicates the poten- 
tial value for naval operational planning. Combinations of 
different conditions of wind, wave, and bottom bathyr.etry can 
be effectively simulated and serve as the basis for determin- 
ing an amphibious landing. 

The open water model is used to forecast sea heights over 
continental shelf and wide stretches of shallow water coast. 

The calculations of wave energy are the result of wind forcing, 
fetch distance, and bottom friction. The model predicts the 
spectral distribution of wave energy. 

The surf model can predict wave shoaling, coastal cur- 
rents and multiple breaker lines. In the model, the wave 
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heights and directions are calculated for a section of beach. 
The breaker heights, breaker types, and location of severe 
breakers are empirically determined based on the wave height 
calculations. This provides a map of the surf zone in two 
dimensions for Navy operations. 

One of the drawbacks of this model is that the model is 
not programmed in an efficient manner. The wavenumber at a 
local water depth is calculated over a hundred times in nearly 
every subroutine for different frequency components. To allow 
for various options chosen and setting up the connectivity 
between subprograms, the integrated model was too large and 
not very well structured. Moreover, the unsteady approach 
to the steady state solution is time consuming and needs 
special care to avoid numerical instability. Using a steady 
state approach and reforming the program structure would 
facilitate the application of this model and improve the 
operational efficiency. Further modeling efforts should be 
directed toward this direction. 

A problem in running the model is the choice of the 
energy bandwidth. The open water model uses energy spectra 
to describe the wave field. The selection of energy bands 
and bandwidths requires knowledge of wave spectra and gravity 
wave theory. The alternative to user choice is to have the 
model calculate energy contribution over the entire wave 
spectra. This would ensure all energy contributions across 
the spectra are considered. However, CPU time would increase 
to more than 6 hours if only 10 energy bands were used. (The 
Navy Spectral Ocean Wave Model uses 15 energy bands.) 
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B. CONCLUSIONS 



This sea, surf and swell model at present is not suitable 
for operational use. The separate regimes, open water and 
nearshore current and surf, all have problems. 

The open water part of the model requires sophisticated 
spectral energy input. The average Aerographer Mate does not 
have the knowledge to properly select input parameters. Only 
when the sea is dominated by a single narrow band swell with 
light wind is the wave energy input fairly easy Cinput wave 
height, direction and period). In addition, the open water 
model was shown to build up the sea height too fast for wind 
speeds of 30 kts raising doubt concerning the wind energy 
transfer part of the model. 

The nearshore current part of the model is very slow to 
converge requiring more than 23 hours to complete the neces- 
sary 1000 iterations. The numerical scheme is therefore 
unsatisfactory for the HP-9845B computer. Recommendations 
include: 1) using a more stable scheme to increase the time 

step and reduce the number of iterations; 2) incorporate an 
implicit numerical scheme to calculate a direct solution of 
the currents such as in Wu and Liu (1985); or 3) use a simple 
one-dimensional current model implying a uniform beach 
alongshore . 

The surf part of the model is probably the most "user 
friendly" and has the most potential. The surf part of the 
model is fast running and gives good solutions for grid point 
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wave heights. However, the width of the surf zone calculation 
is very poorly predicted. This error in turn causes the 
breaker line position, height and angle to be incorrect. 

The incorrect breaker predictions cause the effective surf 
calculation also to be invalid. 
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